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Abstract 

We investigate the long time stability of the Sun-Jupiter-Saturn-Uranus system 
by considering the planar, secular model. Our method may be considered as an 
extension of Lagrange's theory for the secular motions. Indeed, concerning the 
planetary orbital revolutions, we improve the classical circular approximation by 
replacing it with a torus which is invariant up to order two in the masses; therefore, 
we investigate the stability of the elliptic equilibrium point of the secular system for 
small values of the eccentricities. For the initial data corresponding to a real set of 
astronomical observations, we find an estimated stability time of 10 years, which 
is not extremely smaller than the lifetime of the Solar System (~ 5 Gyr). 



1 Introduction 

In this paper we revisit the problem of the stability of the Solar System, at least consid- 
ering (some of) the major planets, in the light of both Kolmogorov's and Nekhoroshev's 
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theories. Some aspects are also related to the theory of Lagrange and Laplace on the 
secular motions of perihelia and nodes of the planetary orbits. 

One of our main aims is to point out the major dynamical and computational diffi- 
culties that arise in the application of Kolmogorov's theorem. In view of this, we attempt 
to apply the Nekhoroshev's theory by trying essentially an extension of Lagrange's theory. 
Although the final results appear to be interesting, our conclusion will be that further 
and more refined investigations are needed. We consider indeed the present paper as the 
beginning of a more comprehensive study of systems with more than two planets in the 
framework of perturbation methods related to the theories above. 

In 1954 Kolmogorov announced his celebrated theorem on the persistence under small 
perturbation of quasi periodic motions on invariant tori of an integrable Hamiltonian 
systems (see [22]). The relevance of that result for the problem of stability of the Solar 
System was pointed out by Kolmogorov himself, and later emphasized in the subsequent 
papers of Moser (see [3JJ]) and Arnold (see pQ). The three papers mentioned above marked 
the beginning of the so called KAM theory. 

However the actual applicability of Kolmogorov's theorem to the planetary system 
encounters two major difficulties, namely: (i) the degeneracy of the Keplerian motion, 
and (ii) the extremely restrictive assumptions on the smallness of the perturbation. 

The former difficulty is related to the elliptic form of the Keplerian orbits. Indeed a 
system including a central body (a star) and n > 1 planets, after elimination of the known 
first integrals, has 3n — 2 degrees of freedom, while only n actions appear in the Keplerian 
part of the Hamiltonian. The way out proposed by Arnold, and inspired by the approach 
of Lagrange and Laplace, was to introduce in the proof two separate time scales for the 
orbital motion and for the secular evolution of the perihelia and of the nodes (see [2] 
and its recent extension in 0). Such an approach has been successfully extended to the 
n + 1-body planetary systems thanks to the work done by Herman and Fejoz (see [TT]). 

Attacking the second difficulty (i.e., the unrealistic requirements on the smallness 
of masses, eccentricities and inclinations of the planets) with purely analytical methods 
seems to be unrealistic. However, some positive results could be attained using computer 
algebra. This means that we explicitly perform a few perturbation steps, thus getting an 
approximation of the wanted invariant torus which is good enough to allow us to apply 
analytical methods. Such an approach (also implementing interval arithmetic) allowed 
some authors to rigorously prove the existence of KAM tori for some interesting problems 
in Celestial Mechanics (see, e.g., [6], [7], [8], [35] and [12]). However, all these works 
consider models having just two degrees of freedom. This because increasing the num- 
ber of independent variables makes the explicit calculation of perturbation steps a big 
challenge, due to the dramatic increase of the number of coefficients to be calculated, 
so that a sufficiently good initial approximation of an invariant torus is hardly obtained. 
For what concerns problems with more than two degrees of freedom, in a few cases only 
the availability of an algorithmic version of Kolmogorov's theorem (see [3J, [TT] and [T8] ) 
allowed us to obtain a good approximation of the invariant tori, although this approach 
is not yet sufficient for a fully rigorous application of the theory. For instance, the con- 
structed solution on a KAM torus has been successfully compared with the real motion 
of the Sun- Jupiter-Saturn system, which can be represented by a model with 4 degrees 
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of freedom (see [37] for all details). Moreover our recent work focuses on a first study of 
the long time stability in a neighborhood of such KAM torus (see [19J). 

Besides the technical difficulty, the results of the numerical explorations have raised 
some doubts concerning the applicability of Kolmogorov's theory to the major bodies 
of our planetary system, namely the Sun and the so called Jovian planets, i.e. Jupiter, 
Saturn, Uranus and Neptune, hereafter we will refer to this model as the SJSUN problem. 
Indeed, the motion of such planetary subsystem has been shown to be chaotic by Sussman 
and Wisdom (see [19]). Murray and Holman provided such an enlightening explanation 
of this phenomenon, that we think it is helpful to briefly summarize some of their results 
as follows (see jlQ] for completeness). 

(a) The chaoticity of the Jovian planets appears to be due to the overlap of some 
resonances involving three or four bodies. An example is given by the resonances 

3ni - 5n 2 - 7n 3 + [(3 - j)g x + Qg 2 + jg 3 ) , with j = 0, 1, 2, 3 , 

where ni stands for the mean motion frequency of the i-th planet, means the (secu- 
lar) frequency of its perihelion argument and the indexes 1, 2, 3 refer to Jupiter, Sat- 
urn and Uranus, respectively. In fact, during the planetary motion each angle cor- 
responding to the resonances above moves from libration to rotation and viceversa. 
Many other resonances analogous to the previous ones are located in the vicinity of 
the real orbit of the SJSUN system, some of them involving also Neptune and the 
frequencies related to the longitudes of the nodes. 

(b) The time needed by these resonances to eject Uranus from the Solar System is 
roughly evaluated to be about 10 18 years. 

(c) By moving the initial semi-major axis of Uranus in the range 19.18-19.35 AU one 
observes some regions that look filled by quasi-periodic ordered motions and other 
regions that are weakly chaotic, i.e., with a Lyapunov time ranging between 2 x 10 5 
and 10 8 years. All the main resonances acting in this region involve the linear 
combination 3n\ — 5n 2 — 7n 3 among the mean motion frequencies of Jupiter, Saturn 
and Uranus. 

(d) The result (c) qualitatively persists also for the planar SJSUN system or when the 
influence of Neptune is neglected. 

(e) Conversely, no chaotic motions are detected in the planar system including the Sun, 
Jupiter, Saturn and Uranus (hereafter, SJSU for shortness) for the same initial 
values of the semi-major axis of Uranus considered at point (c). This suggests that 
the resonances described at point (a) affect observable regions only when combined 
with some effects induced by Neptune or by the mutual inclinations. 

By the way, we note that the resonances involving the linear combination 3n± — 5n 2 — 
7n 3 are clearly related to the approximate ratio 5 : 2 and 7 : 1 between the orbital 
motion of Jupiter and Saturn and of Jupiter and Uranus, respectively. Similarly, the 
ratio 2 : 1 between Uranus and Neptune appears also to be relevant (historically, this 
helped Le Verrier to predict the existence and the location of Neptune). The low order 
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of the latter resonance may explain why the influence of Neptune induces some chaotic 
behaviour, as pointed out in (d) and (e) above. 

The weak chaos in the motion of the Jovian planets makes somehow hopeless the task of 
describing their long-term evolution by a quasi-periodic approximation, as it is provided 
by the KAM theory. Therefore it appears to be more natural to look for exponential 
stability as assured by Nekhoroshev's theory (see [H] and [12] )• Indeed the theorem of 
Nekhoroshev applies to an open set of initial conditions, and states that the stability 
time increases exponentially with the inverse of the perturbation parameter. Our aim is 
to investigate whether the SJSU system may remain close to its current conditions for a 
time that exceeds the lifetime of the system itself; e.g., in our case the age of the Universe, 
which is estimated to be ~ 1.4 x 10 10 years, could be enough. We stress that the rather 
long time reported in (d) concerning the possible dissolution of the SJSUN system seems 
to support our hope. The approach based on Nekhoroshev's theory has been applied 
during the last decades to the case of the Trojan asteroids, producing realistic results 
(see, e.g., [21], [17], [10] and [33]). Concerning the SJSUN system, we expect that a 
combination of both the KAM and the Nekhoroshev theory could prove that the motion 
remains close to an invariant torus for very long times (see [SB] and [19J). 

In the present paper, we restrict our attention to the SJSU planar system, due to the 
huge computational difficulties one encounters during the expansion of the Hamiltonian. 
Indeed, a rather long preliminary work is necessary in order to give the Hamiltonian a 
convenient form for starting more standard perturbation methods (see [35], [36] and [37J). 
We devote sects. [2] and 13.11 to this part of the problem. 

Furthermore, in the line of Lagrange's theory, we focus only on the secular part of 
the Hamiltonian, which is derived in subsect. 13.21 Let us emphasize that all along both 
sects. [2] and [3] we pay a special attention to include all the relevant terms related to the 
three-body mean motion quasi-resonance 3ni — 5n 2 — 7n 3 in view of the remarks reported 
at points (a) and (c) above). 

The secular system turns out to have the form of a perturbed system of harmonic 
oscillators. It can be remarked that the reduction to the plane model makes it possible 
to investigate the stability of the equilibrium using the theorem of Dirichlet. However, 
this is enough if we restrict our attention to the planar secular model, but does not apply 
neither to the spatial case, due to the opposite signs in the secular frequencies of the 
perihelia and of the nodes, nor in the full non secular model. Thus, we think that it is 
also useful to proceed by investigating the stability in Nekhoroshev's sense, since this may 
give indication about the possibility of extending our calculation to more refined models. 
This part is worked out in sect. HI 

Finally, sect, \5\is devoted to the conclusions. 

2 Classical expansion of the planar planetary Hamilto- 
nian 

Let us consider four point bodies Po, P±, P 2 , P3, with masses m , mi, m 2 , m 3 , mutually 
interacting according to Newton's gravitational law. Hereafter the indexes 0, 1, 2, 3 will 



On the stability of the secular . . . planar Sun-Jupiter-Saturn-Uranus system 



5 



correspond to Sun, Jupiter, Saturn and Uranus, respectively. 

Let us now recall how the classical Poincare variables can be introduced so to perform 
a first expansion of the Hamiltonian around circular orbits, i.e., having zero eccentricity. 
We basically follow the formalism introduced by Poincare (see [13] and jS] ; for a modern 
exposition, see, e.g., (30] and [32J ) . We remove the motion of the center of mass by 

using heliocentric coordinates r,- =PqPj , with j = 1, 2, 3 . Denoting by r,- the momenta 
conjugated to r,-, the Hamiltonian of the system has 6 degrees of freedom, and reads 

F(tL) = T {0) (r) + U {0 \r) + T^(r) + U {1) (r) , (1) 



where 



3 
I 

3 

I 

The plane set of Poincare's canonical variables is introduced as 



m 2 m 3 



m mj , , 

Aj = \/y(rn + mj)cij , A-,- = Mj + Uj , 



(2) 



for j = 1,2,3, where aj , ej , Mj and Uj are the semi-major axis, the eccentricity, 
the mean anomaly and the perihelion argument, respectively, of the j-th planet. One 
immediately sees that both and rjj are of the same order of magnitude as the eccentricity 



e 3 . 



Using Poincare's variables (121) , the Hamiltonian F can be rearranged so that one has 

F(A, A, |, V) = (A) + fiF^ (A, A, C, v) , (3) 

where F^ = T^ ' /iF^ 1 ) = TW+f/' 1 '. Here, the small dimensionless parameter /x = 
max{m! / m , m 2 / m , m 3 / m } has been introduced in order to highlight the different 
size of the terms appearing in the Hamiltonian. Let us remark that the time derivative 
of each coordinate is O(ia) but in the case of the angles A. Therefore, according to the 
common language in Celestial Mechanics, in the following we will refer to A and to their 
conjugate actions A as the fast variables, while (£,77) will be called secular variables. 

We proceed now by expanding the Hamiltonian in order to construct the first basic 
approximation of Kolmogorov's normal form. We pick a value A* for the fast actions and 
perform a translation 7a* defined as 

L j = A j -A* j , forj = l,2,3. (4) 

This is a canonical transformation that leaves the coordinates A, £ and n unchanged. 
The transformed Hamiltonian = F o 7a* can be expanded in power series of L, £, 77 
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Table 1: Masses rrij and initial conditions for Jupiter, Saturn and Uranus in our planar 
model. We adopt the UA as unit of length, the year as time unit and set the gravitational 
constant Q = 1 . With these units, the solar mass is equal to (27r) 2 . The initial conditions 
are expressed by the usual heliocentric planar orbital elements: the semi-major axis cij , 
the mean anomaly Mj , the eccentricity Cj and the perihelion longitude Uj . The data are 
taken by JPL at the Julian Date 2440400.5 . 





Jupiter (J — 1) 


Saturn (J = 2) 


Uranus (j = 3) 


m j 

CLj 
Mj 

e 3 


(2tt) 2 /1047.355 

5.20463727204700266 

3.04525729444853654 

0.04785365972484999 

0.24927354029554571 


(2tt) 2 /3498.5 

9.54108529142232165 

5.32199311882584869 

0.05460848595674678 

1.61225062288036902 


(2vr) 2 /22902.98 

19.2231635458410572 

0.19431922829271914 

0.04858667407651962 

2.99374344439246487 



around the origin. Thus, forgetting an unessential constant we rearrange the Hamiltonian 
of the system as 

oo oo oo 

H<T> (L, \,i, R )=n*.L + J2 hfj ] GQ+^££ (L, A, £, rj) , (5) 

ji=2 ji=0j 2 =0 

where the functions j 2 are homogeneous polynomials of degree ji in the actions L_ 
and of degree j% in the secular variables (£, rj) . The coefficients of such homogeneous 

polynomials do depend analytically and periodically on the angles A. The terms 
of the Keplerian part are homogeneous polynomials of degree j\ in the actions L, the 
explicit expression of which can be determined in a straightforward manner. In the latter 
equation the term which is both linear in the actions and independent of all the other 
canonical variables (i.e., n* ■ L) has been separated in view of its relevance in perturbation 
theory, as it will be discussed in the next section. We also expand the coefficients of the 
power series hr?j in Fourier series of the angles A. The expansion of the Hamiltonian 
is a traditional procedure in Celestial Mechanics. We work out these expansions for the 
case of the planar SJSU system using a specially devised algebraic manipulation. The 
calculation is based on the approach described in sect. 2.1 of [35J, which in turn uses the 
scheme sketched in sect. 3.3 of [4"6] . 

The reduction to the planar case is performed as follows. We pick from Table IV of jlH] 
the initial conditions of the planets in terms of heliocentric positions and velocities at the 
Julian Date 2440400.5 . Next, we calculate the corresponding orbital elements with respect 
to the invariant plane (that is perpendicular to the total angular momentum). Finally we 
include the longitudes of the nodes Qj (which are meaningless in the planar case) in the 
corresponding perihelion longitude Uj and we eliminate the inclinations by setting them 
equal to zero. The remaining initial values of the orbital elements are reported in Table [TJ 

Having determined the initial conditions we come to determining the average values 
3) of the semi-major axes during the evolution. To this end we perform a 
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long-term numerical integration of Newton's equations starting from the initial conditions 
related to the data reported in Table dJ After having computed (a* , > a ;s) > we determine 
the values A* via the first equation in fl2]). This allows us to perform the expansion (j3J) of 
the Hamiltonian as a function of the canonical coordinates (L, A, £, rj). In our calculations 
we truncate the expansion as follows, (a) The Keplerian part is expanded up to the 
quadratic terms. The terms hj j 2 include: (bl) the linear terms in the actions L, (b2) all 
terms up to degree 18 in the secular variables (£, 77) , (b3) all terms up to the trigonometric 
degree 16 with respect to the angles A . Our choice of the limits will be motivated in the 
next section. 



3 The secular model 

We look now for a good description of the secular dynamics. A straightforward method 
would be to include in the unperturbed Hamiltonian also the average of the perturbation 
over the fast angles. However, it has been remarked by Robutel (see [13]) that the 
frequencies of the quasi-periodic flow given by this secular Hamiltonian (often called 
of order one in the masses) are quite different from the true ones. The reason lies in the 
effect of the mean motion quasi-resonance 5:2. Therefore we look for an approximation 
of the secular Hamiltonian up to order two in the masses (see, e.g., [29], [31], [15], [35] 
and |34J). To this end wc follow the approach in [371. carrying out two "Kolmogorov-like" 
normalization steps in order to eliminate the main perturbation terms depending on the 
fast angles A . We concentrate our attention on the resonant angles 2Ai — 5A 2 , Ai — 7A 3 
and 3Ai — 5A 2 — 7A3, which are the most relevant ones for the dynamics. Our aim is to 
replace the orbit with zero eccentricity with a quasi periodic one that takes into account 
the effect of such resonances up to the second order in the masses. The procedure is a little 
cumbersome, and requires two main steps that we describe in the next two subsections. 

3.1 Partial reduction of the perturbation 

We emphasize that the Fourier expansion of the Hamiltonian (J3J) is generated just by 
terms due to two-body interactions, and so harmonics including more than two fast angles 
cannot appear. Thus, at first order in the masses only harmonics with the resonant angles 
2Ai— 5A2 and Ai— 7A3 do occur. Actually, harmonics with the resonant angle 3Ai— 5A 2 — 7A3 
are generated by the first Kolmogorov-like transformation, but are of second order in the 
masses, and shall be removed by the second Kolmogorov-like transformation described in 
the next section. 

Let us go into details. We denote by [/] AvK - the Fourier expansion of a function / 
truncated so as to include only its harmonics k- A satisfying the restriction < \k\ < Kp . 
We also denote by (-)\ the average with respect to the angles Ai , A 2 , A3 . The canonical 
transformations are using the Lie series algorithm (see, e.g., [14]). 

We set Kp = 8 and transform the Hamiltonian ([5]) as 'w 02 ^ = exp£ ( 02 ) with 
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the generating function ^Xi (Aj£j?7) determined by solving the equation 

3 „ (02) 6 

E^^-+EK2l A:8 ae>2) = o- (6) 

i*=1 3 



32=0 



Notice that, by definition, ^ ["/] ) = , which assures that equation ([6]) can be solved 
provided the frequencies (n* , n* 2 , n^) are not resonant up to order 8 , as it actually occurs 
in our planar model of the SJSU system. 

The Hamiltonian ~H^° 2 ^ has the same form of in (j3J), with the functions 

A (02) 

replaced by new ones, that we denote by hj generated by the expanding the Lie series 
exp C (02) "HP^> and by gathering all the terms having the same degree both in the fast 
actions and in the secular variables. 

Now we perform a second canonical transformation H^ 02 ^ = exp£^(c>2) T-L^° 2 ^ , where 

the generating function /i X2° 2 \Lli A, Cj r ?) (which is linear with respect to L) is determined 
by solving the equation 

3 o (02) 6 

E» 3 - £ i- + E[C , ] A ,(i.A,f, 2 ) = o. P) 

j=l j j 2 =0 

Again, the Hamiltonian H^ 02 ^ can be written in a form similar to ([3]), namely 

oo oo oo 

?^ 2 > (L, A, £, rj) = n* ■ L + £ ft£*> (L) + /i £ E fc S2 & A, £, 2 ; /*) • (8) 

ii=2 ji=0j2=0 

where the new functions h^f are calculated as previously explained for h^f . Moreover, 
they still have the same dependence on their arguments as h^j 2 in (jSJ). 

If terms of second order in /i are neglected, then the Hamiltonian T-L^° 2 ^ possesses the 
secular 3-dimensional invariant torus L = and £ = 77 = 0. Thus, in a small neighborhood 
of the origin of the fast actions and for small eccentricities the solutions of the system 
with Hamiltonian 'HS 02 ' ) differ from those of its average (% { ° 2) )a by a quantity 0(/i 2 ). 
In this sense the average of the Hamiltonian (JH]) approximates the real dynamics of the 
secular variables up to order two in the masses, and due to the choice Kp = 8 takes into 
account the resonances 5 : 2 between Jupiter and Saturn and 7 : 1 between Jupiter and 
Uranus. 

In this part of the calculation we produce a truncated series which is represented as a 
sum of monomials 

%k,r,A 1 Li 2 Li 3 ^e 2 % 3 vl 1 V?V? c S c. n s(Mi + k 2 \ 2 + Ms) • 

The truncated expansion of W° 2 ^ contains 94 109 751 such monomials. We truncate our 
expansion at degree 16 in the fast angles A and at degree 18 in the slow variables £, 77 (we 
shall justify this choice at the end of the next section). 



On the stability of the secular . . . planar Sun-Jupiter-Saturn-Uranus system 



9 



3.2 Second approximation and reduction to the secular Hamilto- 
nian 

The huge number of coefficients determined till now does not allow us to continue by 
keeping all of them. Therefore, in view that we plan to consider the secular system, we 
perform a partial average by keeping only the main terms that contain the resonant angle 
3Ai — 5A2 — 7A3. More precisely, we first consider the reduced Hamiltonian 

00 

(«< O2 »L. ) A = / i E(C ) K,2;rt> A , 0) 

32=0 

namely we set L = , which results in replacing the orbit having zero eccentricity with 
a close invariant torus of the unperturbed Hamiltonian, and average H^ 02 ^ by removing 
all the Fourier harmonics depending on the angles. Next, we select in H^ 02 ^ the Four- 
ier harmonics that contain the wanted resonant angle 3Ai — 5A2 — 7A3 and add them 
to the Hamiltonian OH]). Finally, we perform on the resulting Hamiltonian the second 
Kolmogorov-like step. With more detail, this is the procedure, which is an adaptation of 
a scheme already used in [35] . 

For {31,32) £ N 2 we select the resonant terms 

AS^A.^ = fi ( hf2 exp [ - i(3A a - 5A 2 - 7A 3 )] exp [i(3A x - 5A 2 - 7A 3 )] + 

H < hff 2 exp [i(3A a - 5A 2 - 7A 3 )] >, exp [ - i(3Ai - 5A 2 - 7A 3 )] . 

(10) 

Actually, this means that in our expression we just remove all monomials but the ones 
containing the wanted resonant angle. Using the selected terms we determine a generating 
function l^Xi^ihiiiV) by solving the equation 

E^^§- + E«(A>e^) = o. (ii) 

3=1 j 32=0 

Here we make the calculation faster by keeping only terms up to degree 9 in (£, 77) , this al- 
lows us to keep the more relevant resonant contributions. Then, still following the proced- 
ure outlined in [33] , we calculate only the interesting part of the transformed Hamiltonian 
exp C 2 (re S ) exp C 2 ( res ) Ti^ 02 ^ , namely we keep in the transformation only the part which 

X2 A* Xi 

is independent of all the fast variables (L, A) . This produces the secular Hamiltonian 
%( sec ) ; which satisfies the formal equation (exp£ 2 (res) exp£ 2 (res) r H ( -° 2 ^) . = "H*- sec ^ + 

* M X2 M Xi ' A 

0{\\L\\) + o(/i 4 ), where 

00 / -, 

^ (sec) (i,!Z) = ME (O, + W \ {xt\c^h^} Lx + 

h= ° ~'~ (12) 

I h=0 J t \ I 32=0 ) ,„ I a 
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Here, we denoted by {•, -} L x and {•, the terms of the Poisson bracket involving only 
the derivatives with respect the variables (L, A) and (£, 77), respectively. 

The Hamiltonian so constructed is the secular one, describing the slow motion of 
eccentricities and perihelia. In view of D'Alembert rules (see, e.g., [B]), it contains only 
terms of even degree and so the lowest order significant term has degree 2. We have 
determined the power series expansion of the Hamiltonian up to degree 18 in the slow 
variables. In order to allow a comparison with other expansions, we reported our results 
up to degree 4 in (£, 77) in appendix [S] 

We close this section with a few remarks which justify our choice of the truncation or- 
ders. The limits on the expansions in the fast actions L have been illustrated at points (a) 
and (bl) at the end of section [21 and they are the smallest ones that are required in order 
to make the Kolmogorov-like normalization procedure significant. Since we want to keep 
the resonant angles 2Ai — 5A 2 , Ai — 7A 3 and 3Ai — 5A 5 — 7A 3 , we set the truncation 
order for Fourier series to 16, which is enough. The choice to truncate the expansion 
at degree 18 in the secular variables (£, rj) is somehow subtler. In view of D'Alembert 
rules the harmonics 2Ai — 5A 2 and Ai — 7A 3 have coefficients of degree at least 3 and 6, 
respectively, in the secular variables. Furthermore, the resonant angle 3Ai — 5A5 — 7A 3 
does not appear initially in the Hamiltonian, but is generated by Poisson bracket between 
the harmonics 2Xi — 5A 2 and Ai — 7A 3 , which produces monomials of degree 9 in (£,77). 

Therefore, we decided to calculate the generating functions xf 2 ^ an d U P to degree 
9 (recall equations ([6]) and ([7])). Finally, in the second Kolmogorov-like step we want to 
keep the secular terms generated by the harmonic 3Ai — 5A5 — 7A 3 , which are produced 
by Poisson bracket between monomials containing precisely this harmonic, and then the 
result has maximum degree 18 in (£,77). This explains the final truncation order for the 
slow variables. 



4 Stability of the secular Hamiltonian model 

The lowest order approximation of the secular Hamiltonian "H^ sec \ namely its quadratic 
term, is essentially the one considered in the theory first developed by Lagrange (see [23] ) 
and furtherly improved by Laplace (see [26], [27] and [28J) and by Lagrange himself 
(see [21], [22]). In modern language, we say that the origin of the reduced phase space 
(i.e., (£, 77) = (0, 0) ) is an elliptic equilibrium point (for a review using a modern formalism, 
see sect. 3 of [1], where a planar model of our Solar System is considered). 

It is well known that (under mild assumptions on the quadratic part of the Hamiltonian 
which are satisfied in our case) one can find a linear canonical transformation (£, 77) = 
V(x,y) which diagonalizes the quadratic part of the Hamiltonian, so that we may write 
n {scc) in the new coordinates as 

On, y) = J2 y (*? + V]) + H 2 ] V) + Hf ] (x, y) + (x, y) + . . . , (13) 
3=0 

where Uj are the secular frequencies in the small oscillations limit and H^J is a homo- 
geneous polynomial of degree 2s + 2 in (x, y) . The calculated values of (cji, U2, w 3 ) in our 
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Table 2: Angular velocities u and initial conditions (x(0),y(0)) for our planar secular 
model about the motions of Jupiter, Saturn and Uranus. The frequency vector u refer to 
the harmonic oscillators approximation of the Hamiltonian (written in (flBl) ) and its 
values are given in r ad /year . 





J = l 


J = 2 


J = 3 


Xj(0) 

%(o) 


-1.1212724892 x 10" 4 
1.5407573458 x 10" 2 
-2.5320810665 x 10" 2 


-1.9688444678 x 10" 5 
-3.0574059274 x 10" 2 
-5.2728862107 x 10" 3 


-1.1134564418 x 10~ 5 
1.1186486403 x 10~ 2 
6.0669645406 x 10~ 3 



case are reported in Table [2J 

Thus, we are led to study the stability of the equilibrium for the Hamiltonian (fT3"|) . As 
remarked in the introduction, perpetual stability in a neighbourhood of the equilibrium 
is assured in our case by Dirichlet's theorem because all frequencies have the same sign, 
that is negative in our case. Actually, a very rough evaluation of the size of the stability 
neighbourhood gives a value about 0.6 times the distance (from the origin) of the actual 
initial data of the planets. Such an estimate should certainly be improved by a more 
accurate calculation, i.e., by determining the stationary points of a function in 6 variables. 
However, we emphasize that our model is just a planar approximation of the true problem. 
If, for instance, one considers the spatial secular problem then the secular frequencies of 
the nodes have a positive sign, so that Dirichlet's theory does not apply any more. Thus, 
we think it is more interesting to investigate the stability of the equilibrium into the light 
of Nekhoroshev's theory. 

4.1 Birkhoff's normal form 

Following a quite standard procedure we proceed to construct the Birkhoff's normal form 
for the Hamiltonian ( fl3|) (see [5]; for an application of Nekhoroshev's theory see, e.g., [.13]). 
This is a well known matter, thus we limit our exposition to a short sketch adapted to 
the present context. 

The aim is to give the Hamiltonian the normal form at order r 

H^(x,y) = Z (§L) + ... + Z r (§L) + T&fey) + J*%(x>y) + • • • , (14) 

where 

*j = \(xj + Vj) for.; = 1,2, 3 (15) 

are the actions of the system, and Z s for s = 0, . . . , r is a homogeneous polynomial of 
degree s/2 + 1 in $ and in particular it is zero for odd s. The un-normalized reminder 
terms Ts , where s > r , are homogeneous polynomials of degree s + 2 in (x, y) . 

We proceed by induction. Assume that the Hamiltonian is in normal form up to a 
given order r, which is trivially true for r = , and determine a generating function 
and the normal form term Z r+ i , by solving the equation 

{x (r+1) , w + rt r Ux,y) = Z r+1 (§L) . (16) 
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Using the algorithm of Lie series transform, we can write the new Hamiltonian as if ( r+1 ) = 
exp £ x (r+i) H^ T \ It is not difficult to show that if( r+1 ) has a form analogous to that written 

in (CHI) with new functions J r J r+1 ' ) of degree s + 2 (where s > r + 1) and the normal form 
part ending with Z r+ \ , which is equal to zero if r is even (see, e.g., [IS]). As usual when 
using the Lie series methods, we denote by (x, y) the new coordinates, so that the normal 
form possesses the approximate first integrals $ given by (fl5|) . By the way, the 
algorithm can be iterated up to the step r provided that the non-resonance condition 

k-u^O Vk G Z 3 such that < \k\ < r + 2 (17) 

is fulfilled. 



4.2 Study of the stability time 

It is well known that the Birkhoff's normal form at any finite order r is convergent in 
some neighbourhood of the origin, but the analyticity radius shrinks to zero when r — > oo . 
Thus, the best one can do is to look for stability for a finite but long time. We use the 
algorithm reported in [20], that we describe here. 

Let us pick three positive numbers R\, R2, R3 and consider a polydisk A e R with center 
at the origin of M 6 defined as 

A e n = {(x,y) G M 6 : x) + y) < g 2 R 2 , j = 1, 2, 3} , 

q > being a parameter. Let Qo = q/2 , and let (x , y ) G A Qo r be the initial point of an 
orbit, so that one has 3>j(0) = (x 2 + y 2 )/2 < g\R 2 j2. Therefore, there is T(g ) > such 
that for |t| < T(q ) we have $(t) < g 2 R 2 /2 , and so also (x(t),y(t)) G A gR . We call T(g ) 
the estimated stability time, and our aim is to give a good estimate of it. 
The key remark is that one has 

00 

*i = & , H^} = {*i > ^ (r) } * > for j = 1, 2, 3 , (18) 

s=r+l 

which holds true for an arbitrary normalization order r. This means that the time deriv- 
ative of $(£) is small, being 0(g r+3 ), so that the time T(qq) may grow very large. The 
basis of Nekhoroshev's theory is that one can choose an optimal value of r as a function 
of £>o letting it to get larger and larger when g — > 0, so that T(g ) grows faster than 
any power of 1/qq . Here we give this argument an algorithmic form, thus producing an 
explicit estimate of T(g ) . 

Let us write a homogeneous polynomial f(x,y) of degree s as 

f(x_,y)= te¥ ' 

|j|+|fc|=s 

where the multiindex notation xpy- = x^xfxfy^y^y^ 3 has been used. We define the 
quantity \f\ R as 

v . ^ . ^ . / j^k^ 

\f\n= 2^ l/i^l-Rf 1 - R 2 2 3 ®ii,fci%!,A:2®i3,fc3 > = y 7 ■ , 7 • (19) 

|j|+|fc]=* V u J 
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We claim that for g > one has 

sup \f(x,y)\<g s \f\ R . (20) 

(x,y)GA eM 

The estimate is checked as follows. In the plane Xi, %ji consider a disk with radius R^ 
Then inside the disk the inequality \xf < Rf + 'O^a- holds true. In fact, after having 
set Xi = Ri cos 9 , yt = Ri sin 9 , one can easily check that | cos^ 0sin 4 0| < 0^,*- . It is 
then straightforward to verify that for a monomial xhj- of degree s one has 

sup \xiy k \ < g s R^Ri^Rf +ks Q nM Q nM Q J3M . 

(x,y)eA sR 

The wanted inequality is just the sum of the contributions of all monomials. 
Using f[2"Uj) and (fTSl) we can estimate 

sup |$,-(a:,y)| < C7 e p+3 |{%,^J 1 }| fl (21) 
(5>»)eA e R - 

for j = 1, 2, 3 and with some C > 1 . In fact, after having set g smaller than the 
convergence radius of the remainder series J-"i r ^ (where s > r), the above inequality is 
true for some C . In our calculation we set C = 2 . 

We come now to the calculation of the estimated stability time. Since $j = g 2 R^/2 , 
we have $j = and, in view of inequality ( l2~Tj) . also 

Thus a majorant of the function g(t) is given by the solution of the equation g = 
B r jg r+2 / Rj . Setting £> as the initial value we conclude that g(t) < 2g for all \t\ < 
T~(go,r), where 

R 2 f 2eo da / 1 \ R 2 
rigo, r) = min — — / = min 1 pr . (22) 

The latter estimate holds true for arbitrary normalization order r. Therefore we select an 
optimal order r opt (^»o) by looking for the maximum over r of r(go, r), thus getting 

T(q ) = maxr(g ,2g ,r) . (23) 

r 

This is the best estimate of the stability time given by our algorithm. 
4.3 Application to the SJSU system 

We apply the algorithm of the previous section to the secular Hamiltonian %( scc ) by ex- 
plicitly performing the construction of Birkhoff 's normal form up to order 30. Meanwhile 
also the first term of the remainder has been stored, so that the estimate for $ is provided. 
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(b) Estimated stability time 

Figure 1: Optimal normalization order r opt and estimated stability time T(g ) evaluated 
according to the algorithm of sect. 14.21 The time unit is the year. See text for more 
details. 
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The calculation of the estimated stability time is performed by setting 



Ri = 2.5558203988 x KT 2 , R 2 = 3.0601862602 x 10~ 2 , R 3 = 1.1223294461 x 10~ 2 



initial data reported in table [2j so that the initial point is on the border of the polydisk 
A gE with g = 1 . 

Finally we proceed to calculating the optimal normalization order r opt (g ) and the 
estimated stability time T(qq) as functions of go in an interval such that the optimal 
normalization order produced by our algorithm is less than 30. The results are reported 
in fig. [U The fast increase of the time when g decreases is evident from the graph. We 
also remark that for g = 1, which corresponds to the initial data for the planets, the 
normalization order is already r opt = 16. This shows that the mechanism of long time 
stability is already active. The estimated time with our algorithm is about 10 7 years for 
go = 1. This seems to be quite short both with respect to the age of the Solar System 
(which is estimated to be ~ 5 x 10 9 years) and with respect to the numerical indications 
(10 18 years). We shall comment on this point in the next section. 

5 Conclusions and outlooks 

In the framework provided by the Nekhoroshev's theorem, the present paper describes the 
first attempt to study the stability of a realistic model with more than two planets of our 
Solar System. As remarked at the end of the previous section we are not yet able to prove 
the stability for a time comparable to the age of our planetary system, even restricting 
ourselves to consider just the secular part of a planar approximation including the Sun, 
Jupiter, Saturn and Uranus. Nevertheless, we think that our results is meaningful in that 
it indicates that the phenomenon of exponential stability in Nekhoroshev sense may play 
an effective role for the Solar System, at least for the biggest planets. On the other hand, 
we stress that our result is not dramatically far from the goal of proving stability for the 
age of the Solar System: such a time is reached for g Q ~ 0.7 . By the way, it may be worth 
to note that a similar result, with the same value of the radius, has been found in [2D] 
where the spatial problem for the Sun- Jupiter-Saturn system is considered. Such a value 
of go appears to be not so small, especially if one recalls the rough estimates based on 
the first purely analytical proofs of the KAM theorem: in order to apply them to some 
model of our planetary system, the Jupiter mass should be smaller than that of a proton. 
Improvements are surely possible, and the relatively short history of the applications 
of the Nekhoroshev's type estimates to Celestial Mechanics has shown definitely more 
remarkable improvements than the one required here (e.g., compare [15] with [21]). 

Some drawbacks are immediately evident. The most relevant one is that the estim- 
ate in ( I22p actually assumes that the perturbation constantly forces the worst possible 
evolution. This is clearly pessimistic, and justifies the striking difference with respect to 
the indication given by the numerical integrations. On the other hand, general perturb- 
ation method are essentially based on estimates that are often very crude. The explicit 



(24) 



These values have been calculated 
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calculation of normal forms and related quantities allows us to significantly improve our 
results, but the price is either a bigger and bigger computer power or more and more 
refined methods. 

The natural question is whether there is a way to improve the present result. Our 
approach suggests that a better approximation of the true orbit could help a lot. This 
can be obtained, e.g., by first establishing the existence of a KAM torus close to the 
initial conditions of the planets, and then proving the stability in Nekhoroshev's sense 
in a neighbourhood of the torus that contains the initial data. Such an approach has 
been attempted in [T9] for the Sun-Jupiter-Saturn case considering the full system, i.e., 
avoiding the approximation of the secular model. In that case the number of coefficients to 
be handled is so huge that the calculation can actually be performed only by introducing 
strong truncations on the expansions; this might artificially improve the results. Thus, 
some new idea is necessary, and this will be work for the future. 
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A Expansion of the secular Hamiltonian of the planar 
SJSU system up to order 2 in the masses and 4 in 
eccentricities 

Our secular model is represented by the Hamiltonian 7i^ sec \ which is defined in ( I12p . Here, 
we limit ourselves to report the expansion of H (scc) up to de gree 4 in (£, rj) . Therefore, as 
a consequence of the D'Alembert rules, the terms related to the quasi-resonance 3Ai — 
5A 2 — 7A 3 do not give any contribution to the coefficients listed below. Thus, the following 
expansion of the rhs of f[T2l actually takes into account just /^{h^^x + }\ ( reca U 

that %( sec ) contains just terms of even degree in its variables (£, rf) ). The calculation of 
the functions an d ft-cu is performed how it has been explained in subsect. 13.11 
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